function [Q,errorMes,model,setupFilter,outFilter] = objectFuncQML(paramsInputValues,setupFilter)
Q         = NaN;
model     = [];
outFilter = [];
% Check if the lower and upper bounds are violated
if isfield(setupFilter,'transformParamsOn') == 1
    paramsInputValues = unTransformParams(paramsInputValues,...
        setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues);
else
    if any(paramsInputValues < setupFilter.lowerBoundsValues) ...
            || any(paramsInputValues > setupFilter.upperBoundsValues)
        errorMes = 1;
        return;
    end
end

% We turn paramsInputValues into a struct 
for i=1:length(setupFilter.selectParams)
    name = setupFilter.selectParams(i,1);
    paramsInput.(name{1}) = paramsInputValues(i,1);
end

% Accounting for calibrated coefficients, we unold paramsInput by to get params 
params = paramsInput;
% The setdiff implies that we do not overwrite elements in params based on calibrateparams
namesCalibrate = setdiff(fieldnames(setupFilter.calibrateParams),setupFilter.selectParams);
for i=1:length(namesCalibrate)
   name = namesCalibrate(i);  
   params.(name{1}) = setupFilter.calibrateParams.(name{1}); 
end

% Setting the size of the measurement errors
params.stdMhours =  params.stdMhours*params.ScalingMacro;
params.stdMwage  =  params.stdMwage*params.ScalingMacro;
params.stdMdc    =  params.stdMdc*params.ScalingMacro;
params.stdMinf   =  params.stdMinf*params.ScalingMacro/2;

params.stdMr1    =  params.stdMr1*params.ScalingYield;
params.stdMr4    =  params.stdMr4*params.ScalingYield;
params.stdMr12   =  params.stdMr12*params.ScalingYield;
params.stdMr20   =  params.stdMr20*params.ScalingYield;
params.stdMr28   =  params.stdMr28*params.ScalingYield;
params.stdMr40   =  params.stdMr40*params.ScalingYield;
params.stdSurvey =  params.stdSurvey*params.ScalingYield;

if ~isfield(params,'ZI')
    params.ZI = (1-params.THETA+params.ETA*params.THETA)*(params.ETA-1)*params.XI/...
        ((1-params.XI)*(1-params.THETA)*(1-params.XI*params.BETTA*params.MUZss^(1-params.CHI*(1-params.CHI0)-params.CHI0)));
end

if ~isfield(params,'ALFA') && (isfield(params,'U0d') || isfield(params,'U0'))
    params.ALFA = getALFA(params);
    if abs(params.ALFA) > 500
        errorMes = 1;
        return;
    end
elseif ~isfield(params,'U0d') && isfield(params,'U0')
    params.U0d = getU0d(params);
elseif isfield(params,'U0d') && ~isfield(params,'U0')
    params.U0 = getU0d(params);
end

%% Solve the model
[model,errorMes] = perturbationDSGEmodel(params,setupFilter.appOrder,...
    setupFilter.setupModel,setupFilter.maturities,setupFilter.MEXon);
if errorMes == 1
    disp('DSGE model: No unique and stable solution')
    Q = NaN;
    errorMes = 1;
    return
end
model.inclSurveys = setupFilter.inclSurveys;
model.pruningOn   = setupFilter.pruningOn;
model.heteroShocks = zeros(1,model.ne);
%% Running the filter
% Setting the co-variance matrix for the measurement errors
if setupFilter.inclSurveys == 1
    dimy  = setupFilter.numMacro+length(model.maturities)+length(model.matConShortRate); 
else
    dimy  = setupFilter.numMacro+length(model.maturities);
end
Rw        = zeros(dimy,dimy);
Rw(1,1)   =  params.stdMhours^2;
Rw(2,2)   =  params.stdMwage^2;
Rw(3,3)   =  params.stdMdc^2;
Rw(4,4)   =  params.stdMinf^2;
Rw(5,5)   =  params.stdMr1^2;
Rw(6,6)   =  params.stdMr4^2;
Rw(7,7)   =  params.stdMr12^2;
Rw(8,8)   =  params.stdMr20^2;
Rw(9,9)   =  params.stdMr28^2;
Rw(10,10) =  params.stdMr40^2;
if setupFilter.inclSurveys == 1
    Rw(11:dimy,11:dimy) = params.stdSurvey^2*eye(dimy-10);
end
model.Rw= Rw;
varX     = reshape((eye(model.nx^2)-kron(model.hx,model.hx))\reshape((model.eta*model.eta'),model.nx^2,1),model.nx,model.nx);
vec_xfxf = reshape(varX,model.nx^2,1);
if setupFilter.pruningOn == 1
    % Initializing the filter
    % We use the pruned perturbation approximation
    if model.appOrder == 1
        xf0       = zeros(model.nx,1);
        xAll0     = xf0;
        PxAll0    = ones(model.nx+0*model.nx1,model.nx+0*model.nx1)*1D-5;
        PxAll0(1:model.nx,1:model.nx) = reshape(vec_xfxf,model.nx,model.nx);
    elseif model.appOrder == 2
        xf0       = zeros(model.nx,1);
        xs0       = (eye(model.nx)-model.hx)\(model.Hxx*vec_xfxf +0.5*model.hss);
        xAll0     = [xf0;xs0(1:model.nx1,1)];
        PxAll0    = ones(model.nx+1*model.nx1,model.nx+1*model.nx1)*1D-5;
        PxAll0(1:model.nx,1:model.nx) = reshape(vec_xfxf,model.nx,model.nx);
    elseif model.appOrder == 3
        xf0       = zeros(model.nx,1);
        xs0       = (eye(model.nx)-model.hx)\(model.Hxx*vec_xfxf +0.5*model.hss);
        xrd0      = xf0*0;
        xAll0     = [xf0;xs0(1:model.nx1,1);xrd0(1:model.nx1,1)];
        PxAll0    = ones(model.nx+2*model.nx1,model.nx+2*model.nx1)*1D-5;
        PxAll0(1:model.nx,1:model.nx) = reshape(vec_xfxf,model.nx,model.nx);
    end
    % Running the filter
    [outFilter,errorMes]=CDKF_dsge(model,setupFilter,setupFilter.data(:,1:dimy),xAll0,PxAll0,Rw);
else    
    % Initializing the filter
    x0      = (eye(model.nx)-model.hx)\(model.Hxx*vec_xfxf +0.5*model.hss);
    P0      = varX;
    
    % Running the filter
    [outFilter,errorMes]=CDKF_dsge(model,setupFilter,setupFilter.data(:,1:dimy),x0,P0,Rw);
end

if errorMes ~= 0
    return;
end
if setupFilter.SMMwithCSreg == 2
    % Adding term premia for nominal yields to model
    TP_version = 2;
    model = getTermPremia_loadings(model,setupFilter,TP_version);
    termPrem = getTermPremia(model,outFilter.xhat);
    % Risk-adjusted CS regression loadings
    % Convert into forwards
    fwdPrem = nan(size(termPrem.TP));
    for i = 1:size(termPrem.TP,1)-1
        fwdPrem(i,:) = (i+1)*termPrem.TP(i+1,:) - i*termPrem.TP(i,:);
    end
    % TPxhat and hence fwdPrem are scaled by 40000 to get risk premia in
    % annualized basis points in getTermPremia.
    % Annualized yields in are not scaled. Hence, we need to undo the scaling of TP
    % to get basis points, but we still need to preserve the annualization i.e. the multiplication
    % by 4.
    numBoots = 0;
    CS_riskAdj = CampbellSchillerRegression_riskAdjust(setupFilter.yieldsAll,...
        termPrem.TP'/40000*4,fwdPrem'/40000*4,setupFilter.CSregress_m,numBoots);
end
%% Unconditional Moments for Shrinkage 
%addMerror = 0;
if setupFilter.numSim > setupFilter.burnIn
    % Moments computed by simulation
    outSim      = simPerturbation(model,setupFilter.shocks(:,1:setupFilter.numSim));
    if setupFilter.pruningOn == 1
        %yieldsAllSim   = max(getYields(model,outSim.xAll_cu),0);
        yieldsAllSim   = getYields(model,outSim.xAll_cu);
    else
        yieldsAllSim   = getYields(model,outSim.x_cu);
    end
    yieldsAllSim         = 4*yieldsAllSim;

    % With measurement noise
    yieldsAllSimWithNoise = 4*yieldsAllSim + model.params.stdMr1*setupFilter.noiseYieldsAll(1:size(yieldsAllSim,1),1:size(yieldsAllSim,2));
    
    pos_l                = find(strcmp({'$l_t$'},model.labely));
    pos_w                = find(strcmp({'$w_t$'},model.labely));
    %pos_c_ba1           = find(strcmp({'$c_{t-1}$'},model.labelx));
    %pos_myz             = find(strcmp({'$\mu _{z,t}$'},model.labelx));
    %pos_c               = find(strcmp({'$c_t$'},model.labely));
    pos_dc               = find(strcmp({'$dc_t$'},model.labely));
    pos_pai              = find(strcmp({'$\pi_t$'},model.labely));
    ySimAsInDataWithNoise= zeros(setupFilter.numMacro+length(model.maturities),setupFilter.numSim);
    
    ySimAsInDataWithNoise(1,:)    = outSim.y_cu(pos_l,:)-model.g0(pos_l,1) + model.params.stdMhours*setupFilter.noise(1,1:size(outSim.y_cu(pos_l,:),2));
    ySimAsInDataWithNoise(2,:)    = outSim.y_cu(pos_w,:)-model.g0(pos_w,1) + model.params.stdMwage*setupFilter.noise(2,1:size(outSim.y_cu(pos_w,:),2));
    %ySimAsInDataWithNoise(3,:)   = 4*(outSim.y_cu(pos_c,:)-(outSim.x_cu(pos_c_ba1,:)+model.g0(pos_c_ba1,1))+outSim.x_cu(pos_myz,:)+log(model.params.MUZss)) ...
    %                               + model.params.stdMdc*setupFilter.noise(3,1:size(outSim.y_cu(pos_c,:),2));
    ySimAsInDataWithNoise(3,:)    = 4*outSim.y_cu(pos_dc,:) + model.params.stdMdc*setupFilter.noise(3,1:size(outSim.y_cu(3,:),2));   %gives position 3, which is ok.
    ySimAsInDataWithNoise(4,:)    = 4*outSim.y_cu(pos_pai,:)+ model.params.stdMinf*setupFilter.noise(4,1:size(outSim.y_cu(pos_pai,:),2));
    ySimAsInDataWithNoise(setupFilter.numMacro+1:setupFilter.numMacro+length(model.maturities),:) ...
                                  = yieldsAllSimWithNoise(model.maturities,:);
    
    % Without noise added
    ySimAsInData                  = zeros(setupFilter.numMacro+length(model.maturities),setupFilter.numSim);
    ySimAsInData(1,:)             = outSim.y_cu(pos_l,:)-model.g0(pos_l,1);
    ySimAsInData(2,:)             = outSim.y_cu(pos_w,:)-model.g0(pos_w,1);
    %ySimAsInData(3,:)            = 4*(outSim.y_cu(pos_c,:)-(outSim.x_cu(pos_c_ba1,:)+model.g0(pos_c_ba1,1))+outSim.x_cu(pos_myz,:)+log(model.params.MUZss));
    ySimAsInData(3,:)             = 4*outSim.y_cu(pos_dc,:);
    ySimAsInData(4,:)             = 4*outSim.y_cu(pos_pai,:);
    ySimAsInData(setupFilter.numMacro+1:setupFilter.numMacro+length(model.maturities),:) ...
                                  = yieldsAllSim(model.maturities,:);                 
                     
    burnIn                           = setupFilter.burnIn;
    outFilter.yieldsAllSim           = yieldsAllSim;
    outFilter.yieldsAllSimWithNoise  = yieldsAllSimWithNoise;
    outFilter.outSim                 = outSim;
    outFilter.ySimAsInData           = ySimAsInData;
    outFilter.ySimAsInDataWithNoise  = ySimAsInDataWithNoise;
    if any(any(abs(outSim.y_cu')> 500)) || any(any(isnan(outSim.y_cu')))
        disp('DSGE model: explodes in simulation')
        errorMes = 1;
        return
    else
        if setupFilter.SMMwithCSreg == 2
            [modelSMM,CSmodel] = momentsSMM(ySimAsInData(:,burnIn+1:end)',...
                yieldsAllSim(:,burnIn+1:end)',setupFilter.nyMom,setupFilter.CSregress_m,1);
            % We add the covariances for the risk-adjusted CS moments to the model moments 
            select   = setupFilter.maturites_CS/4;
            T        = size(CS_riskAdj.Ydata,1);
            tmp_yx   = (CS_riskAdj.Ydata(1:T,select)-repmat(mean(CS_riskAdj.Ydata(1:T,select),1),T,1)).*CS_riskAdj.Xdata(1:T,select);
            % For checking
            %tmp_xx   = (CS_riskAdj.Xdata(1:T,select)-repmat(mean(CS_riskAdj.Xdata(1:T,select),1),T,1)).*CS_riskAdj.Xdata(1:T,select);
            %[CS_riskAdj.CSbetta(2,setupFilter.maturites_CS/4) mean(tmp_yx,1)./mean(tmp_xx,1)]
            modelSMM.dataMoments = [modelSMM.dataMoments mean(tmp_yx,1)];
        else
            [modelSMM,CSmodel] = momentsSMM(ySimAsInData(:,burnIn+1:end)',...
                yieldsAllSim(:,burnIn+1:end)',setupFilter.nyMom,setupFilter.CSregress_m,setupFilter.SMMwithCSreg);
        end
        outFilter.modelMoments = modelSMM.dataMoments;
        
        % Extra moments to evaluate the fit of the model
        outFilter.resMoments.CS_loadings  = CSmodel.CSbetta(2,:);
        outFilter.resMoments.meanYields   = mean(yieldsAllSim(:,burnIn+1:end),2);
        outFilter.resMoments.stdYields    = std(yieldsAllSim(:,burnIn+1:end),0,2);
        outFilter.resMoments.meanForEstim = mean(ySimAsInData(1:setupFilter.nyMom,burnIn+1:end),2);
        outFilter.resMoments.stdForEstim  = std(ySimAsInData(1:setupFilter.nyMom,burnIn+1:end),0,2);
    end
elseif setupFilter.tau == 0
    % Moments computed in closed form based on the pruned approximation
    outFilter.resMoments   = momentsClosedForm(model,setupFilter);
    %scaling                = outFilter.resMoments.scaling;
    if setupFilter.SMMwithCSreg == 2
        % We add the covariances for the risk-adjusted CS moments to the model moments 
        select   = setupFilter.maturites_CS/4;
        T        = size(CS_riskAdj.Ydata,1);
        tmp_yx   = (CS_riskAdj.Ydata(1:T,select)-repmat(mean(CS_riskAdj.Ydata(1:T,select),1),T,1)).*CS_riskAdj.Xdata(1:T,select);
        outFilter.resMoments.modelMoments = [outFilter.resMoments.modelMoments mean(tmp_yx,1)];
    end
    outFilter.modelMoments = outFilter.resMoments.modelMoments;    

end
   
%% The objective value
T                      = setupFilter.T;
outFilter.momCondition = outFilter.modelMoments-setupFilter.dataMoments;
outFilter.SMMpart      = setupFilter.lambdaSMM*outFilter.momCondition*setupFilter.W*outFilter.momCondition';
outFilter.QMLpart      = sum(outFilter.LogL(setupFilter.numInitial+1:end))/(T-setupFilter.numInitial); 
outFilter.LogL         = outFilter.LogL; 
Q                      = -(outFilter.QMLpart - outFilter.SMMpart);
outFilter.Q            = Q;
% if setupFilter.numSim > setupFilter.burnIn
%     outFilter.outSim              = outSim;
%     outFilter.ySimAsInData        = ySimAsInData;
%     outFilter.ySimAsInDataNoNoise = ySimAsInDataNoNoise;
% end
  
    
end

